/*
 * Gaussian_Reuse_Parallel.c
 *
 *  Created on: Apr 29, 2009
 *      Author: yzhang8
 */

#include "Gaussian.h"
#include "Gaussian_Reuse_Parallel.h"

inline int min(int a, int b);
void* compute_row_barrier(void* s);

int volatile col;
pthread_barrier_t thread_barrier, total_barrier;
pthread_mutex_t mutex;
pthread_cond_t cond;

int gaussian_elimination_reuse_parallel() {
	time_t time_sec_start, time_sec_finish;
	time_sec_start = time(NULL);
	pthread_t* p_threads = malloc(thread_num * sizeof(pthread_t));
	pthread_attr_t attr;
	pthread_attr_init(&attr);
	pthread_cond_init(&cond, NULL);
	if (pthread_barrier_init(&thread_barrier, NULL,thread_num) != 0) {//for the threads
		//printf("failed to init barrier");
		return EXIT_FAILURE;
	}
	if (pthread_barrier_init(&total_barrier, NULL,thread_num + 1) != 0) {//include itself
		//printf("failed to init barrier");
		return EXIT_FAILURE;
	}
	if (pthread_mutex_init(&mutex, NULL) != 0) {
		//printf("failed to get mutex");
		return EXIT_FAILURE;
	}
	thread_row_range* thread_arg = malloc(thread_num * sizeof(thread_row_range));
	int i, j, thread_first_row, thread_row_size, useful_thread_num;
	double *row_swap = malloc(size * sizeof(double)), b_swap;
	char buffer[40];
	for (col = 0; col < size - 1; col++) {// go through each column
		pthread_mutex_lock(&mutex);//the sequential step
		//printf("col %d: main thread gets lock\n", col);
		//fflush(stdout);
		int pivot_row = col;
		for (i = col + 1; i < size; i++) {//find the largest
			if (matrix_A[pivot_row][col] < matrix_A[i][col]) {
				pivot_row = i;
			}
		}
		if (pivot_row != col) {//exchange rows
			b_swap = vector_B[col];
			vector_B[col] = vector_B[pivot_row];
			vector_B[pivot_row] = b_swap;
			for (j = col; j < size; j++) {
				row_swap[j] = matrix_A[col][j];
				matrix_A[col][j] = matrix_A[pivot_row][j];
				matrix_A[pivot_row][j] = row_swap[j];
			}
		}
		useful_thread_num = (size - col - 1) / thread_row_threshold + 1;//calculate how many threads should be used
		if (useful_thread_num > thread_num) {
			useful_thread_num = thread_num;
		}
		thread_first_row = col + 1;
		thread_row_size = (size - col - 1) / useful_thread_num;
		for (i = 0; i < useful_thread_num; i++) {
			thread_arg[i].first_row = thread_first_row;
			thread_arg[i].thread_id = i;
			if (i != useful_thread_num - 1) {//not the last thread
				thread_arg[i].last_row = thread_first_row + thread_row_size - 1;
				thread_first_row = thread_arg[i].last_row + 1;
			} else {
				thread_arg[i].last_row = size - 1;
			}
		}
		for (; i < thread_num; i++) {//for the unused thread
			thread_arg[i].thread_id = -i; //set it as useless
		}
		pthread_mutex_unlock(&mutex);
		//printf("col %d: main thread release lock\n", col);
		//fflush(stdout);
		if (col == 0) {//create all threads
			for (i = 0; i < thread_num; i++) {//run threads
				pthread_create(&p_threads[i], &attr, compute_row_barrier, (void*) &thread_arg[i]);
			}
		} else {
			pthread_mutex_lock(&mutex);
			//printf("col %d: main thread gets lock again\n", col);
			//fflush(stdout);
			pthread_cond_signal(&cond);
			//printf("col %d: main thread unblock condition\n", col);
			//fflush(stdout);
			pthread_mutex_unlock(&mutex);
			//printf("col %d: main thread release lock again\n", col);
			//fflush(stdout);
		}
		pthread_barrier_wait(&total_barrier);//wait for the threads to finish
		//printf("col %d: main thread passed total barrier\n", col);
		//fflush(stdout);
		sprintf(buffer, "\nGaussian after proccess col %d\n ", col);
		print_Gaussian(buffer);
		//fflush(stdout);
	}
	pthread_mutex_lock(&mutex);
	//printf("col %d: main thread gets lock again\n", col);
	//fflush(stdout);
	pthread_cond_signal(&cond);
	//fflush(stdout);
	//printf("col %d: main thread unblock condition after finish\n", col);
	pthread_mutex_unlock(&mutex);
	//printf("col %d: main thread release lock again\n", col);
	//fflush(stdout);
	for (i = 0; i < useful_thread_num; i++) {//wait for all thread to complete calculate the rows
		pthread_join(p_threads[i], NULL);
	}
	pthread_barrier_destroy(&thread_barrier);
	pthread_barrier_destroy(&total_barrier);
	pthread_cond_destroy(&cond);
	pthread_mutex_destroy(&mutex);
	time_sec_finish = time(NULL);
	printf("\nGaussian elimination used %d sec\n", (int) (time_sec_finish - time_sec_start));
	free(p_threads);
	free(thread_arg);
	return EXIT_SUCCESS;
}

void* compute_row_barrier(void* s) {
	thread_row_range* range = (thread_row_range*) s;
	int row, j, jj, jsize, rsize, rr;
	double* pivot_row = malloc(size * sizeof(double));
	double* row_ratio = malloc(size * sizeof(double));
	while (1) {//all threads need to exit after the barrier
		if (range->thread_id == 0) {//thread leader
			pthread_mutex_lock(&mutex);//the sequential step
			//printf("col %d: thread 0 gets lock\n", col);
			//fflush(stdout);
		}
		//printf("col %d: thread %d wait for thread barrier\n", col, range->thread_id);
		//fflush(stdout);
		pthread_barrier_wait(&thread_barrier);
		//printf("col %d: thread %d past thread barrier\n", col, range->thread_id);
		//fflush(stdout);
		if (col == size - 1) {
			break;
		}
		if (range->thread_id >= 0) {//do
			for (j = col; j < size; j++) {
				pivot_row[j] = matrix_A[col][j];
			}
			for (row = range->first_row; row <= range->last_row; row++) {
				row_ratio[row] = matrix_A[row][col] / pivot_row[col];
			}
			for (row = range->first_row; row <= range->last_row; row += 50) {
				rsize = min(row + 50, range->last_row + 1);
				for (j = col; j < size; j += 50) {
					jsize = min(j + 50, size);
					for (rr = row; rr < rsize; rr++) {
						for (jj = j; jj < jsize; jj++) {
							matrix_A[rr][jj] -= pivot_row[jj] * row_ratio[rr];
						}
					}
				}
			}
			for (row = range->first_row; row <= range->last_row; row++) {
				vector_B[row] -= vector_B[col] * row_ratio[row];
			}
		}
		//printf("col %d: thread %d wait for total barrier\n", col, range->thread_id);
		//fflush(stdout);
		pthread_barrier_wait(&total_barrier);
		//printf("col %d: thread %d past total barrier\n", col, range->thread_id);
		//fflush(stdout);
		if (range->thread_id == 0) {//thread leader
			//printf("col %d: thread 0 gets to wait \n", col);
			//fflush(stdout);
			pthread_cond_wait(&cond, &mutex);
			//printf("col %d: thread 0 wake up from wait \n", col);
			//fflush(stdout);
			pthread_mutex_unlock(&mutex);//
			//printf("col %d: thread 0 release lock after wait\n", col);
			//fflush(stdout);
		}
	}
	//printf("thread %d  exits\n", range->thread_id);
	//fflush(stdout);
	free(pivot_row);
	free(row_ratio);
	pthread_exit(0);
}

inline int min(int a, int b) {
	return ((a > b) ? b : a);
}
